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The  Use  of  Matrix  Displacement  Method  for 
Vibrational  Analysis  of  Structures 


by 

William  0.  Hughes1 


Department  of  Mechanical  Engineering  and 
Robotics  Institute 
Carnegie-Mellon  University 
Pittsburgh,  PA  15213 
May,  1981 

ABSTRACT  A  study  of  the  matrix  displacement  method  for  modeling  the  vibrations  of 
structures  is  presented  in  this  report.  The  model  can  analyze  both  the  free  and 
forced  vibrations  of  a  structure.  Static  loading  on  a  structure  is  treated  as  a 
special  case  of  the  forced  vibration  analysis. 
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1  Introduction 

A  study  of  the  matrix  displacement  method  for  modeling  the  vibrations  of  structures  is  presented  in  this 
report  The  model  can  analyze  both  the  free  and  forced  vibrations  of  a  structure.  Static  loading  on  a  structure 
is  treated  as  a  special  case  of  the  forced  vibration  analysis. 


A  brief  review  of  the  Finite  Element  Method  and  its  present  use  is  first  given.  This  is  followed  by  a 
discussion  of  the  methodology  of  the  matrix  displacement  approach  and  a  description  of  the  specific  model 
used.  Examples  of  the  use  of  the  model  to  analyze  the  frequencies  and  mode  shapes  of  the  free  and  forced 
response  of  a  beam  structure  and  the  static  deflections  of  a  beam  structure  are  shown  and  compared  with  the 
closed  form  solutions.  Finally,  ways  of  extending  the  model  to  a  more  complicated  structure,  a  turbine  blade, 
are  discussed.  Conclusions  are  then  drawn.  ('-  - 

2  The  Finite  Element  Method  ••  Fundamental  Concepts  and  Applications 

There  are  many  methods  available  today  which  perform  the  analysis  of  structures.  For  example,  in  one 
method  the  structure  is  described  by  differential  equations.  The  differential  equations  are  then  solved  by 
analytical  or  numerical  methods.  Another  method  of  analysis  is  the  finite  element  method  (FEM). 

In  this  method,  the  structure  is  idealized  into  an  assembly  of  discrete  structural  elements,  each  having  an 
assumed  form  of  displacement  or  stress  distribution.  The  complete  solution  is  then  obtained  by  assembling 
these  individual,  approximate,  displacement  or  stress  distributions  in  a  way  satisfying  the  force  equilibrium 
equations,  the  constitutive  relationships  of  the  material,  the  displacement  compatibility  between  and  within 
the  elements  and  the  boundary  conditions  of  the  structure. 

Methods  based  on  discrete  element  idealization  have  been  used  extensively  in  structural  analysis.Thc  early 
pioneering  works  of  Turner,  et  al.,  in  1956  [1],  and  Argyris  in  1960  [2]  led  to  the  application  of  this  method  to 
static  and  dynamic  analysis  of  aircraft  structures.  Other  fields  of  structural  engineering,  such  as  nuclear 
reactor  design  and  ship  construction  have  since  employed  this  method. 

Nor  is  the  idea  of  discrete  elements  limited  in  use  to  structural  analysis  only.  The  fundamental  concept  of 
the  finite  element  method  is  that  any  continuous  quantity,  such  as  displacements,  temperature,  or  pressure, 
can  be  approximated  by  a  finite  number  of  elements.  Thus,  this  approach  can  be  used  to  solve  problems  in 
heat  flow,  fluid  dynamics,  electro-magnetics,  fracture  mechanics  and  seepage  flow  to  name  just  a  few  other 
areas  of  usage. 


The  representation  of  a  continuous  structure  by  structural  elements  of  finite  size  results  in  large  systems  of 
algebraic  equations.  A  convenient  way  of  handling  these  sets  of  equations  is  by  the  use  of  matrix  algebra, 
which  also  has  the  advantage  of  being  ideally  suited  for  computations  on  high-speed  digital  computers.  For 
this  reason,  expressions  such  as  "matrix  methods  of  structural  analysis"  are  sometimes  used  to  describe  the 
method.  More  common  though  is  the  term  "finite  element  method",  which  emphasizes  the  discretisation  of 
the  structure. 

The  finite  element  method  actually  encompasses  three  classes  of  matrix  methods  of  structural  analysis.  The 
first  is-  the  displacement  (or  stiffness  method),  where  the  displacements  of  the  nodes  are  considered  the 
unknowns.  The  correct  set  of  displacements  results  from  satisfying  the  equations  of  force  equilibrium.  The 
second  method  is  the  force  (or  flexibility)  method.  Here  the  nodal  forces  are  the  unknowns  and  are  found  by 
satisfying  the  conditions  of  compatible  of  deformations  of  the  members.  The  third  class  of  matrix  method  is 
the  mixed  method,  which  is  a  combined  force-displacement  method. 

One  last  comment  on  the  finite  element  method  in  general  is  necessary.  An  error  is  introduced  into  the 
solution  of  the  original  problem  as  soon  as  the  continuous  structure  is  replaced  by  discrete  elements.  This 
error  remains,  even  when  the  discrete  element  analysis  is  performed  exactly.  In  general  this  error  is  reduced 
by  increasing  the  number  of  discrete  elements,  thereby  decreasing  the  element  size  and  thus  giving  a  better 
idealization  of  the  continuous  structure.  Zienkiewicz,  Brotton  and  Morton  [3]  suggest  that  the  us^r  may 
determine  the  limits  of  his  error  by:  "(a)  comparison  of  finite  element  calculations  with  exact  solutions  for 
cases  similar  to  his  specific  problem;  (b)  a  ’convergence  study’  in  which  two  or  more  solutions  are  obtained 
using  progessively  finer  subdivisions  and  the  results  plotted  to  establish  their  trend  or  (c)  using  experience  of 
previous  calculations  as  a  guide  to  the  treatment  of  the  specific  problem."  Further  information  on  matrix 
structural  analysis  and  the  finite  element  method  may  be  found  in  many  sources.  [4-11] 

3  Explanation  of  the  Model 

The  following  discussion  is  divided  into  three  sections.  Firstly  the  equations  of  motion  will  be  stated. 
Secondly,  the  matrix  displacement  method  for  solving  such  equations  will  be  described.  Finally  some  specific 
aspects  of  the  particular  model  being  used  will  be  discussed. 

3.1  Equations  of  Motion 

The  motion  of  a  vibrating  system,  consisting  of  mass  and  stiffness,  of  n  degrees  of  freedom  can  be 
represented  by  n  differential  equations  of  motion.  These  equations  of  motion  may  be  obtained  by  Newton’s 
second  law  of  motion,  by  ^grange's  equation  or  by  the  Influence  Coefficients  method.  Since  the  equations 
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of  motion,  in  general,  are  not  independent  of  each  other,  a  simultaneous  solution  of  these  equations  is 
required  to  calculate  the  frequencies  of  the  system. 


The  matrix  equation  for  the  free  vibration  case  is: 
[K  -  «2M]  pC]  =  [0] 


(1) 


where 


[K]  represents  the  stiffness  matrix  of  the  structure, 

[M]  represents  the  inertial  (mass)  matrix  of  the  structure, 

u  represents  the  set  of  eigenvalues  of  the  equations 

corresponding  to  the  set  of  natural  frequencies, 

[X]  represents  the  set  of  eigenfunctions  of  the  equations 

corresponding  to  the  set  of  displacements 


For  the  free  vibration  case  the  set  of  forces  is  just  zero. 

The  matrix  [K-w2M]  is  called  the  impedance  matrix. 

The  matrix  equation  for  the  forced  vibration  case  is: 

[K  -  «2  M]  PC]  =  (PJ  (2) 

where  [P]  represents  the  set  of  forces  on  the  structure,  and 

wf  is  the  driving  or  forcing  frequency. 

The  other  terms  are  as  previously  defined. 

Inspection  of  equations  (1)  and  (2)  reveals  that  neither  contain  damping  terms.  This  is  because  structures 
of  immediate  concern  have  very  low  damping  (~1  x  10'4  critical  damping). 

An  excellent  treatment  on  the  dynamics  of  structures  is  Clough  and  Pcnzicn  [14]. 


3.2  The  Matrix  Displacement  Method 

An  outline  of  the  application  of  the  matrix  displacement  method  in  finite  element  analysis  for  the  solution 
of  dynamic  problems  follows.  A  similar  outline  is  given  by  Zienkiewicz,  et  al.  [3]  for  static  analysis. 

1.  Input 

a.  Idealization  of  the  problem 

The  continuous  structure  is  divided  into  a  number  of  elements.  These  elements  are 
connected  at  common  nodal  points  or  nodes.  It  is  at  these  nodes  that  the  value  of  the 
continuous  quantity  (displacement)  is  to  be  determined. 

b.  Preparation  of  the  data  for  the  structure 

The  geometry  of  the  structure  is  defined  by  assigning  coordinates  to  the  nodal  points.  The 
physical  properties  of  the  elements  (dimensions,  material  parameters)  are  inputted. 

c.  Preparation  of  the  load  data 

The  loads  to  be  applied  to  each  element  or  node  are  defined. 

d.  Preparation  of  the  boundary  conditions  or  constraints 

The  prescribed  constraints  on  the  degrees  of  freedom  and  boundary  conditions  are  stated. 

2.  Processing 

a.  Element  Formulation 

The  stiffness  and  inertial  matrices  for  each  element  are  determined  by  the  approximate 
relationships  and  the  corresponding  loads  are  calculated. 

b.  Assembly  of  the  structure 

The  summation  of  the  elemental  matrices  to  form  structural  stiffness,  inertial  and  load 
matrices  is  performed. 

c.  Reduction  of  equations 

The  boundary  conditions  and  constraints  in  terms  of  certain  specified  displacements  are 
introduced,  thereby  reducing  the  number  of  equations  to  be  solved. 

d.  Solution,  of  simultaneous  equations 

The  solution  of  the  eigen  problem  of  equation  (1)  or  (2)  results  in  the  natural  frequencies  of 
the  structure  (eigenvalues)  and  the  modal  shapes  or  displacements  of  the  nodes 
(eigenfunctions). 


e.  Calculation  of  stresses 

If  required,  the  elemental  stresses  could  be  calculated  from  the  nodal  displacements  and 
elemental  stiffness. 

3.  Output 

The  results  of  the  solution  to  the  eigenvalue  problem  and  the  stress  calculation  are  presented  in  an  easily 
interpreted  form. 

3.3  Specific  Aspects  of  Model  • 

This  section  is  concerned  with  specific  aspects  of  the  model.  The  element  and  its  formation  M  be 
discussed  first  Information  concerning  the  computer  code  and  its  subroutines  will  then  be  given. 

1.  Element  Formulation 

The  element  chosen  for  the  model  is  the  beam  element  which  is  given  by  Przemicmiecki  [7].  This  element 
was  chosen  so  as  to  allow  direct  comparison  of  results  with  known  solutions  (see  section  4).  The  beam 
element  is  a  two  node  element  The  model  allows  the  nodes  to  have  either  three  degrees  of  freedom  (x  and  y, 
translational  and  rotation  about  z,  Le.  motion  confmed  to  a  plane)  or  six  degrees  of  freedom  (x,y,z 
translational,  rotation  about  x,y,z,  i.e.  the  general  case). 

Fig.  1  shows  the  beam  element  The  following  forces  act  on  the  beam: 

•  axial  forces  sx  and  &7 

•  shearing  forces  s^  s3,  Sg,  and  s^ 

•  bending  moments  s5,  s6,  sn,  and  s^ 

•  and  twisting  moments  (torques)  s4  and  s10. 

The  location  and  positive  directions  of  these  forces  are  also  given  in  Fig.  1.  The  corresponding 
displacements  Up  U2,. . .  U12  will  be  taken  to  be  positive  in  the  positive  direction  of  these  forces. 

Each  element  has  its  own  set  of  physical  parameters.  For  the  beam  clement  these  parameters  arc:  Young’s 
modules,  cross-sectional  area,  moment  of  inertia  about  the  y  and  z  axis,  Poisson's  ratio,  mass  density,  and 
length  (along  x  axis).  All  of  these  parameters  are  inputted  directly  except  for  the  length  which  is  computed 
from  the  inputted  coordinates  of  the  nodes.  * 
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Figure  1:  The  beam  element  and  its  forces,  after  Przemieniccki  [7] . 

The  model  performs  calculations  for  cither  the  free  or  forced  vibration  case.  To  perform  such  calculations 
requires  the  calculation  of  the  structural  stiffness  and  inertial  matrices,  along  with  information  of  the  loading 
and  boundary  conditions  of  the  structure.  The  effect  of  constraining  a  degree  of  freedom  is  to  strike  out  the 
corresponding  rows  and  columns  of  the  stiffness,  interial  and  load  matrices. 

The  stiffness  matrix  for  a  beam  clement  is  shown  in  Fig.  2.  The  shear  deformation  parameters  «I>y  and 
can  be  taken  as  zero.  This  matrix  may  be  obtained  in  various  ways,  two  of  which  arc  the  influence  coefficients 
method  and  the  variational  method,  which  are  outlined  in  Appendices  I  and  II. 

The  inertial  matrix  for  the  beam  clement  is  shown  in  Fig.3.  This  matrix  is  obtained  by  the  same  methods  as 
the  stiffness  matrix,  as  described  in  Appendices  I  and  II. 

Liepcns  [13]  gives  a  third  way  of  calculating  the  stiffness  and  inertial  matrices. 

The  structural  matrix  for  both  stiffness  and  inertia  is  obtained  by  superposition  of  the  individual  elemental 
matrices.  Actual  superposition  occurs  only  when  degrees  of  freedom  arc  common  to  more  than  one  element 

2.  Computer  Coding 

The  computer  code  itself  contains  ten  subroutines,  called  by  the  main  program,  entitled  VIBRAT.  A  brief 
explanati  on  of  the  subroutines  will  now  be  given. 


INPUT  -  This  subroutine  asks  the  user  tor  the  necessary  information  which  is  needed  to  assemble  the 
structure.  Information  such  as:  free  or  forced  case,  number  of  elements,  coordinates  of 
nodes,  physical  parameters,  structural  loading,  and  constrained  degrees  of  freedom  are 
inputted  in  this  section. 

CONECT  -  This  subroutine  establishes  the  geometry  of  the  model.  It  determines  the  distances  between 
adjacent  nodes  of  the  structure. 

KM  AT  -  This  subroutine  calculates  the  elemental  stiffness  matrix  for  each  element  and  then  assembles  the 
structural  stiffness  matrix  from  them. 

MMAT  -  This  is  similar  to  KMAT  only  here  the  mass  or  inertial  matrices  are  calculated. 

EIGEN  -  This  subroutine  is  called  for  the  free  vibration  case.  The  purpose  of  it  is  to  calculate  the  eigenvalues 
(natural  frequencies)  and  eigenvectors  (mode  shapes)  of  equation  (1).  This  subroutine  calls 
two  other  subroutines:  EIGZF,  an  IMSL  routine  which  actually  does  the  solving,  and 
CLAMPR,  which  determines  which  degrees  of  freedom  are  constrained. 

SOLVE  -  This  subroutine  is  called  for  the  forced  vibration  case.  This  routine  solves  equation  (2)  for  the 
displacement  This  subroutine  also  calls  two  other  subroutines:  LEQT1F,  an  IMSL  routine 
which  does  the  solving,  and  CLAMPR,  which  determines  the  proper  degrees  of  freedom  to 
be  constrained. 

REMARK  -  is  a  subroutine  whose  purpose  is  to  explain  the  use  of  the  main  program  VIBRAT  and  its 
subroutines.  Information  on  the  nomenclature  and  file  structure  used  can  be  found  in 
REMARK.  The  user  of  the  model  is  recommended  to  refer  to  REMARK  if  he  has  any 
questions  on  the  computer  code  used  in  this  model. 


The  code  for  all  of  these  routines  may  be  found  in  Appendix  III. 


4  The  Model:  Examples  and  Accuracy 

This  section  presents  various  examples  of  use  of  the  model.  The  examples  chosen  represent  five  types  of 
possible  problems.  They  are: 

1.  free  vibration  of  a  fixed-free  uniform  beam 

2.  free  vibration  of  a  fixed-fixed  uniform  beam 

3.  forced  vibration  of  a  fixed-free  uniform  beam 

4.  static  deflection  of  a  fixed-free  uniform  beam 

5.  static  deflection  of  a  fixed-free  non-uniform  beam. 


The  accuracy  of  each  example  is  discussed. 
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Figure  2:  Stiffness  Matrix  of  Beam  Element  of  Figure  1  [After  PrzmienieckiJ. 
[The  sheer  deformation  parameters  $y  and  <&  can  be  considered 

to  be  zero.] 


The  first  four  examples  use  the  geometric  and  material  values  listed  in  Table  1. 


Parameter 


Value 


Total  Beam  Length  (L) 

Young's  Modulus  (E) 
Cross-Sectional  Area  (A) 

Moment  of  I nertia  about  Z-Axis  (Iz) 
Moment  of  Inertia  about  Y-Axis  (Iy) 
Poisson’s  Ratio  (?) 

Mass  Density  (p) 


25.0 

27.8  x  106 
2.0 
0.2 
0.7 
0.305 
0.283 


inches 

pounds  force/inches2 

inches2 

inches4 

inches4 

pounds  mass/inches3 


Table  1:  Uniform  Beam  Properties 


Figure  3:  Consistent  Mass  Matrix  for  a  Beam  Element 
(A  fter  Przemieniecki  [7D , 


Example  1:  Free  Vibration  of  a  Fixed-Free  Uniform  Beam 


Figure  4:  Example  1:  Fixed-Free  Uniform  Beam. 

Table  2  summarizes  the  results  for  this  problem,  using  one,  two.  and  five  elements.  It  is  clear  that 


increasing  the  number  of  elements  increases  the  accuracy  of  the  results,  and  this  supports  the  statements  of 
Zienkiewicz  given  earlier. 

The  natural  frequencies  calculated  by  the  model  are  compared  with  the  closed  form  solution  obtained 
from  the  partial  differential  equation  of  the  continuous  system.  For  the  fixed-free  case  the  closed  form 
solutions  are: 


Axial  u>  =  — — — >/-  —  where  n  =  1, 3, 5, . . .  (3) 

2L  p 

where  1  +  cos  aL  cosh  aL  =0 

i=Y  or  Z  (4) 

Torsional  «  =  — — \J—  9.  —  where  n=l,  3,  5, . .  .G= — - —  (5) 

2L  p  2(1  +  p) 

Thus  from  Table  2,  one  can  see  that  by  using  just  five  elements,  the  model  gives  ten  transverse  modes,  two 
axial  modes,  and  two  rotational  modes,  the  frequencies  of  which  are  all  within  5%  of  the  exact  solutions. 
Again,  clearly  greater  accuracy  of  results  and  more  (higher)  modes  may  be  accomplished  by  increasing  the 
number  of  elements. 

Diagrams  of  the  mode  shapes  for  the  first  five  bending  modes  (in  Y)  and  the  first  four  axial  modes  (along 
X)  are  given  in  Figs.  5  and  6.  The  model  shapes  agree  with  the  closed  form  predictions  in  every  case. 

4.2  Example  2:  Free  Vibration  of  a  Fixed-Fixed  Uniform  Beam 

In  this  example  the  beam  is  held  fixed  on  both  ends.  See  Figure  7  .  Table  3  shows  the  calculated  and  exact 
values  for  the  axial  mode  natnral  frequencies.  The  accuracy  is  similar  to  that  of  example  1. 

4.3  Example  3:  Forced  Vibration  of  a  Fixed-Free  Uniform  Beam 

In  this  example  (Figure  8),  the  beam  is  subjected  to  a  harmonically  varying  load  P(t)  of  amplitude  P  and 
circular  frequency,  Figure  9  is  a  plot  of  the  magnitude  in  the  transverse  direction  of  the  free  end  node.  As 
expected,  as  o>f  approaches  a  natural  frequency  (those  found  in  example  1).  a  resonance  condition  occurs 
resulting  in  very  large  magnitudes  of  deflection.  The  expression  for  the  amplitude  of  response  A  is  given  by 
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Axial 

Mode 

Calculated  Natural 
Frequency  (rad/sec) 

Exact  Natural 
Frequency  (rad/sec) 

% 

Error 

1 

24,874 

24,470 

1.7 

2 

52,186 

48,940 

6.6 

3 

83,933 

73,410 

14.3 

4 

117 ; 570 

97,880 

20.1 

Table  3:  Calculated  and  Exact  Natural  Frequencies  in  Axial  Mode. 
Calculated  value  used  five  element  model,  for  Example  2. 


— ^ — D  (6) 

K 

represents  the  static  deflection, 

equals  the  ratio  of  the  forcing  frequency  to  natural  frequency, 
dynamic  magnification  factor  equal  to  l/(l-/82) 

Analysis  of  the  calculated  amplitude  in  terms  of  the  dynamic  magnification  factor  agrees  with  equation  (6) 
in  those  frequency  regions  dominated  by  just  one  natural  frequency. 

4.4  Example  4:  Static  Deflection  of  a  Fixed-Free  Uniform  Beam 

By  letting  the  driving  frequency,  be  zero  in  the  forced  vibration  option,  the  model  is  able  to  solve  static 

deflection  problems.  Figure  11  shows  the  deflection  of  the  beam  under  the  static  loading  of  example  4.  The 

model’s  calculations,  using  just  five  elements  are  within  2%  of  the  exact  beam  theory  results.  The  deflection 

and  slope  at  the  end  of  the  beam  arc  given  by  the  expressions: 

A  =  PL3/3EI 
0  =  PL2/2EI 


A  = - 0-5-= 

K(l-j82) 

where  Po/K 

P 

D 


Figure  8:  Example  3:  Fixed-Free  Uniform  Beam  With  Dynamic  Load 


Values  calculated  using  these  expressions  are  compared  with  the  model  results  in  Table  4 . 

4.5  Examples  5:  Static  Deflection  of  a  Fixed-Free  Non-Uniform  Beam 

Until  now,  all  the  examples  have  dealt  with  uniform  beams.  Example  5  is  an  example  taken  from  Laursen 
[11].  Laursen  solves  the  problem  in  three  differential  ways:  by  the  moment-area  method,  by  the  conjugate 
beam  method,  and  by  Newmark’s  method.  The  solution  for  displacement  and  slope  at  the  free  end  is  given 
as: 

A  =  -0.457  inches 
9  =  -0.0041  radians 


The  model  gives  identical  results. 

A  sketch  of  the  deflection  is  shown  in  Figure  13. 

The  purpose  of  the  previous  five  examples  is  to  illustrate  the  use  and  application  of  the  model  to  a  variety 
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jdc  versus  Forcing 


Frequency  for  Example  3. 


Figure  11:  Static  Deflection  of  a  Uniform  Beam,  Example  4  . 

of  eases.  Other  eases  of  a  more  complicated  nature  could  have  been  solved  as  easily,  however  these  examples 
give  the  user  some  insight  into  the  accuracy  of  the  solution  obtained.  They  also  indicate  that  very  accurate 
results  arc  obtained  by  the  model  with  relatively  few  elements.  In  general,  for  a  more  complicated  structure 
more  elements  will  be  required  to  obtain  an  accurate  model.  Techniques  for  handling  more  complex 
structures  arc  discussed  in  the  next  section. 


Deflection  x  10  ( inches) 


A  (inches) 

0  (radians) 

Exact 

-9.37  x  10"4 

-5.62  x  10'S 

-4 

_e 

Calculated 

-9.50  x  10 

-5.69  x  10  ‘ 

% 

1.4 

1.2 

Tabic  4:  Calculated  and  Exact  Values  of  Deflections  for  Example  4  # 


t 

N 


5  kips 
B 


I  «  S00  in.4 
6  ft 


4 


/  -  500  in.4 
Oft 


Figure  12:  Example  5:  Static  Deflection  of  a  Fixed-Free 
Non-Uniform  Beam . 

[After  Laursen}. 

5  The  Extension  of  the  Model  to  Model  A  Turbine  Blade 

An  example  of  a  more  complicated  structure  which  might  be  of  vibrational  interest  to  an  engineer  is  a 
turbine  blade.  The  equations  of  motion  for  a  beam  in  bending  vibration  is  a  fourth-order  differential 
equation,  whose  solution  is  easily  found.  The  solution  for  a  non-uniform  and  asymmetrical  beam  is  much 
more  complicated.  A  tapered,  pre-twisted  turbine  blade  with  airfoil  cross-section  might  be  modeled  as  such  a 
beam. 


The  differential  equations  for  combined  flapwisc  bending,  chordwisc  bending  and  torsion  of  a  twisted 
non-uniform  blade  arc  derived  by  Houbolt  and  Brooks  [16].  The  solutions  of  these  equations  for  the 
continuous  system  have  not  been  found.  Thus  the  analysis  of  such  structures  arc  limited  to  special  cases 


which  solutions  are  obtainable,  or  to  approximate  solutions.  Various  techniques  of  an  analytical  and  iterative 
nature  such  as  the  Myklestad  method,  Holzcr  method,  Stodala  method,  Raylcigh-Ritz  method,  transmission 
matrix  method,  and  the  Runge-Kutta  method  have  been  studied  [14].  A  few  typical  examples  are  given  in  the 
references  [15,17-20]. 

The  application  of  the  model  presented  in  this  report  to  the  turbine  blade  would  be  a  very  useful  tool  to  the 
engineer  and  his  study  of  the  blade’s  free  and  forced  vibrations. 

The  model  allows  each  element  to  have  its  own  set  of  geometric  and  physical  parameters.  Thus  neither  the 
non-uniformity  or  tapering  of  the  blade  would  lead  to  any  modeling  problems.  However  the  airfoil  shape  of 
the  blade  would  not  have  the  same  torsional  stiffness  as  a  beam.  Thus  the  first  adaptation  to  the  model 
needed  would  be  to  correctly  compute  the  torsional  stiffness  for  an  airfoil  shape  and  input  this  into  the  model 
rather  than  using  that  which  the  model  computes. 

There  is  another  problem  which  arises  from  the  twisting  and  geometry  of  the  turbine  blade.  The  natural 
frequencies  of  such  a  blade  are  coupled  frequencies  with  the  mode  shapes  consisting  in  general  of  transverse 
motion  coupled  with  torsion.  The  coupling  is  dependent  upon  the  degree  of  pre-twist  and  the  ratio  of  depth 
taper  to  width  taper.  For  a  given  blade,  coupling  becomes  stronger  with  increasing  pre-twist  and  with 
increasing  width  to  depth  taper  ratio. 

The  simulation  of  this  coupling  in  the  model  could  be  accomplished  by  either  introducing  it  through  the 
element  itself  or  through  the  geometry  of  the  structure.  The  first  way  implies  changing  the  element  from  a 
beam  element  to  a  new  element  This  new  element  could  be  derived  from  a  variational  method  (see 
Appendix  II)  applied  to  the  differential  equations  for  the  blade  equations  derived  by  Houbolt  and  Brooks 
[16].  The  ideal  of  coupling  through  the  geometry  of  the  structure  implies  the  use  of  additional  beam 
elements.  Part  of  these  elements  would  be  used  to  form  the  center  of  stiffness  for  the  blade  which  would  now 
be  a  curve  rather  than  the  straight  line  used  thus  far.  Other  elements  could  extend  at  right  angles  from  this 
curve.  These  elements  would  act  primarily  as  lumped  masses  and  form  the  curve  representing  the  center  of 
mass  of  the  blade. 

Modeling  a  turbine  blade  with  this  model  would  require  some  additional  work  to  implement  the  ideas 
presented  in  this  section.  However  the  matrix  displacement  method  used  is  a  very  powerful  one  and  the  use 
of  the  model  and  extensions  of  it  are  applicable  to  a  wide  range  of  problems  in  vibrational  analysis  of 
structures.  Building  a  library  of  elements  would  greatly  extend  the  usefulness  of  the  existing  model,  and 
additionally,  the  introduction  of  clement  rotation  would  lead  to  further  improvement 
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6  Conclusion 

This  report  primarily  concerns  itself  with  three  topics: 

1.  the  explanation  of  the  matrix  displacement  method  for  use  in  vibrational  analysis  of  structures, 

2.  specific  examples  showing  the  variety  and  accuracy  of  the  method,  and 

3.  possible  extensions  of  the  model  to  allow  for  application  to  an  even  wider  variety  of  problems. 

The  model  presented  here  currently  allows  for  only  one  type  of  element,  the  beam  element  It  has  been 
shown  that  by  using  just  a  few  beam  elements  very  accurate  results  of  frequencies  and  modal  shape  are 
obtained  for  beam-like  structures.  Creating  a  library  of  element  types  would  allow  the  user  even  greater 
flexibility.  The  accuracy  of  the  model  using  these  new  elements  should  be  comparable  to  that  presented  here. 
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I.  Appendix  I  Influence  Coefficient  Method 

One  method  of  obtaining  the  stiffness  matrix  is  the  influence  coefficient  method.  This  method  is  widely 
used  in  structural  analysis  with’ static  loadings  [10,11].  There  arc  both  stiffness  and  flexibility  influence 
coefficients  :  only  the  stiffness  influence  coefficients  will  be  considered  here. 

The  stiffness  coefficients  for  an  clement  are  found  by  alternatively  constraining  all  degrees  of  freedom  but 
one  and  displacing  this  one  by  a  unit  amount  The  resulting  forces  on  the  other  degrees  of  freedom  are  the 
stiffness  coefficients.  That  is  Kij  is  the  force  or  couple  corresponding  to  degree  of  freedom  t  due  to  the  unit 
displacement  of  degree  of  freedom  j.  In  Fig.  14  a  prismatic  element  of  length  1,  area  A,  moment  of  inertia 
about  the  Z  axis  I,  and  modulus  of  elasticity  E  with  three  degrees  of  freedom  per  node  is  shown. 
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Figure  15:  Stiffness  matrix  of  prismatic  elements  of  Figure  14  , 

Comparison  of  Fig.  2  and  15  shows  that  the  matrix  of  Figure  15  is  contained  within  the  matrix  of  Figure  2. 
In  Fig.  15,  each  node  has  three  degrees  of  freedom,  in  Fig.  2  there  are  six  degrees  of  freedom  per  node. 

The  inertial  (or  mass)  matrix  may  be  calculated  similarly.  The  mass  influence  coefficients  would  represent 
the  mass  inertia  force  acting  at  a  degree  of  freedom  due  to  a  unit  acceleration  of  another  degree  of  freedom. 
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il.  Appendix  li  Variational  Method 

Another  method  of  computing  elemental  stiffness  matrices  is  the  variational  or  energy  method  commonly 
used  in  finite  clement  programs.  The  outline  presented  here  largely  follows  that  of  Gallagher  (8|. 

The  principle  of  minimum  potential  energy  furnishes  a  variational  basis  for  the  formulation  of  the  element 
stiffness  matrix.  The  potential  energy  ( ir p  of  a  structure  is  given  by  the  strain  energy  (U)  plus  the  potential  of 
the  external  work  V  (V  =  -W  ).  The  theorem  of  potential  energy  is:  of  all  displacements,  satisfying  the 
boundary  conditions,  those  that  satisfy  the  equilibrium  conditions  make  the  potential  energy  assume  a 
stationary  (extreme)  value.  Thus 

Wp  =  U  +  V 

8v  =  5U  +  5V  =  0 
P 

And  for  stable  equilibrium,  is  a  minimum. 

8im  =  52U  +  52V  >  0  •  (9) 

P 

The  change  in  strain  energy  density  due  to  the  change  in  strain  caused  by  a  virtual  displacement  (8e)  is  given 
by 

5  (dlD  =  a  8e  (10) 

Where  a  is  the  equilibrium  stress  state  prior  to  the  application  of  the  virtual  displacement.  The  stress-strain 
law  is 

o  =  [E]e  -  [E]e  tat  (11) 

where  [E]  is  called  the  material  stiffness  matrix,  a  matrix  of  elastic  constants.  For  simplicity,  let  ther  e.  be  no 
initial  strain.  Substitution  of  (11)  into  (10)  yields 

8(dlj  =  e  [E]5e  (12) 

Integration  between  zero  and  the  strain  e,  corresponding  to  a,  gives 

dU  =  — ? —  e  (E]e  *  (13) 

2 


(7) 

(8) 


and  integration  over  the  volume  of  the  element  results  in 


U  = 


-/ 


e[E]e  d(vol) 


vol 


The  variation  of  U  is 


-I 


8U  =  J  e[E]  8e  d(vol) 

vol 

The  potential  of  the  applied  loads  is 
V  =  -21FA  -J  T.  Uds 

S 

a 

where  Ft  represents  point  forces,  and  T  are  traction  forces  on  the  surface.  The  variation  of  V  is 

SV  =  *2FtSAt-  J  f  •  SUds 

So 

Using  the  minimum  potential  energy  theorem  (equation  8)  results  in 


/ 


e[E]5e  d(vol)  +  SF 


.“.7  t- 


5  Uds  =  0 


vot 


(14) 


(15) 


(16) 


(17) 


(18) 


In  the  finite  element  matrix,  the  displacements,  [A],  are  written  as  a  polynomial  matrix  times  a  vector  of 
parameters  in  the  assumed  displacement  field. 


[A]  =  [P]  [a] 


(19) 


[P]  evaluated  at  the  node  gives  a  matrix  [B],  consisting  of  constants.  Thus 

[\<xJ  =  lfilN  (20) 

Inverting  to  find  [a]  in  (20)  and  substitution  into  (19)  leads  to 
[A]  =  [PUB*1]  [Anodes] 

=  [N]  Anodes  (21) 


where  N  is  the  shape  function.  The  shape  function  Nt  has  the  quality  that  it  is  equal  to  1  when  evaluated  at 
the  geometric  coordinates  of  the  point  at  which  A(  is  defined  and  is  equal  to  zero  at  all  other  dcgrccs-of- 
freedom  A  ,  J  *t. 


Figure  16:  Axial  element  cross-sectional  area  A,  modulus  E 
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The  result  is  also  contained  in  the  stiffness  matrices  shown  in  Figures  2  and  15. 

The  inertial  (or  mass)  matrix  can  also  be  calculated  by  use  of  this  method.  The  variational  approach  leads 


[N]  [p]  [N]‘dVol 


(28 


where  [p  ]is  the  material  mass  density  matrix.  Since  the  shape  functions  used  here  are  the  same  as  those 
used  for  the  stiffness  calculation  the  result  is  called  the  consistent  mass  matrix.  A  consistent  mass  matrix  is 
more  accurate  than  a  lumped  mass  approach  [12]. 


III.  Appendix  III  Compute 
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